# This is the function that creates the p values for a given set of parameters 
powerTest <- function(screens = 1, sets = 2, items = 6, total.items = 256, 
                      sd = 3.5, respondents = 3000, mean = 5, effect.size = 0.5, 
                      reps = 1000) {
  
  # this function runs one example
  expPValue <- function(screens = 1, sets = 2, items = 6, total.items = 256, 
                        sd = 3.5, respondents = 3000, mean = 5, effect.size = 0.5) {
    
    prob <- (1 / total.items)
    prob.in.set <- (1-(1- prob)^items)
    total.sets <- respondents * sets * screens
    
    x <- rbinom(n = total.sets, prob = prob.in.set, 
                size = 1)
    
    base.utility <- rnorm(mean = mean, sd = sd, n = total.sets)
    
    final.utility <- base.utility 
    final.utility[x==1] <- final.utility[x==1] + effect.size
    
    final.utility[final.utility>10] <- 10
    final.utility[final.utility<0] <- 0
    
    output <- summary(lm(final.utility~ x))
    return(output$coefficients[2, 4])
  }
  p.values <- rep(NA, reps)
  for(ii in 1:length(p.values)) {
    p.values[ii] <- expPValue(screens = screens ,
              sets = sets,
              items = items,
              total.items = total.items,
              sd = sd,
              respondents = respondents,
              mean = mean,
              effect.size = effect.size)  
  }
  return(p.values)
}


# set a seed if you want to be able to replicate the results 
# as there is a random element to these calculations
set.seed(2373247)

final.option.1 <- powerTest(screens = 3, sets = 2, items = 6, total.items = 256, 
                            sd = 3.5, respondents = 5000, mean = 5, effect.size = 0.5, 
                            reps = 1000)

prop.table(table(final.option.1<0.05))


final.option.1.con <- powerTest(screens = 3, sets = 2, items = 6, total.items = 256, 
                                sd = 3.5, respondents = 2000, mean = 5, effect.size = 0.5, 
                                reps = 1000)

prop.table(table(final.option.1.con<0.05))

final.option.1.power.25 <- powerTest(screens = 3, sets = 2, items = 6, total.items = 256, 
                            sd = 3.5, respondents = 5000, mean = 5, effect.size = 0.25, 
                            reps = 1000)

prop.table(table(final.option.1.power.25<0.05))


actual.vals <- powerTest(screens = 3, sets = 2, items = 6, total.items = 257, 
                            sd = 1.31, respondents = 4908, mean = 5.96, effect.size = 0.46, 
                            reps = 10000)

prop.table(table(actual.vals<0.05))

